Limited associations between MHC diversity and reproductive success in a bird species with biparental care

Abstract The selective pressure from pathogens on individuals can have direct consequences on reproduction. Genes from the major histocompatibility complex (MHC) are central to the vertebrate adaptive immune system and pathogen resistance. In species with biparental care, each sex has distinct reproductive roles and levels of investment, and due to a trade‐off with immunity, one can expect different selective regimes acting upon the MHC of each parent. Here, we addressed whether couples combine each other's variation at MHC loci to increase their breeding success. Specifically, we used a 23‐year dataset from a barn owl population (Tyto alba) to understand how MHC class Iα and IIβ functional divergence and supertypes of each parent were associated with clutch size and fledging success. We did not detect associations between MHC diversity and supertypes with the clutch size or with the fledging success. In addition, to understand the relative contribution from the MHC of the genetic parents and the social parents, we analyzed the fledging success using only a cross‐fostered dataset. We found several associations of weak‐to‐moderate effect sizes between the father's MHC and fledging success: (i) lower MHC‐Iα divergence in the genetic father increases fledging success, which might improve paternal care during incubation, and (ii) one and two MHC‐IIβ DAB2 supertypes in the social father decrease and increase, respectively, fledging success, which may affect the paternal care after hatching. Furthermore, fledging success increased when both parents did not carry MHC‐IIβ DAB1 supertype 2, which could suggest conditional effects of this supertype. Although our study relied on a substantial dataset, we showed that the associations between MHC diversity and reproductive success remain scarce and of complex interpretation in the barn owl. Moreover, our results highlighted the need to incorporate more than one proxy of reproductive success and several MHC classes to capture more complex associations.

1 Corresponding author: dianaferreira58@gmail.com 2 Co-last authors  Supplementary material and methods

S1. MHC amplification and genotyping
DNA was extracted from blood with DNeasy Blood and Tissue commercial kit (Qiagen) following the manufacturer's protocol.For the present study, we investigated the exon-3 from the two loci of the MHC-Iα gene and the exon-2 from the two loci of the MHC-IIβ gene (i.e., MHC-IIβ DAB1 and MHC-IIβ DAB2).These regions encode for the highly polymorphic peptide-binding region (PBR) responsible for antigen recognition and are the most studied regions in avian research (Minias et al., 2018).All individuals were co-amplified for both loci of MHC-Iα while loci of MHC-IIβ were amplified independently, both using a dual barcode strategy designed for high-throughput sequencing.Forward and reverse primers were modified by adding a random label composed of three nucleotides (NNN) and an eight-nucleotide barcode sequence to the 5′-end.The three random nucleotides (NNN) were added to increase the diversity because the MiSeq technology requires nucleotide diversity to properly identify clusters on the flow cell (Gaigher et al., 2016).
Briefly, MHC-Iα and MHC-IIβ fragments were amplified in a final volume of 25 μL containing 1 × buffer Gold, 2.0 mM MgCl2 Gold, 1× Q solution (Qiagen), 0.2mM dNTPs, 0.25 μM of each primer, and 0.5 U AmpliTaq Gold (Applied Biosystems).PCR conditions included an initial denaturation step at 95°C for 10 min, 32 cycles of denaturation at 95°C for 30 s, annealing at 68°C (60°C for MHC-IIβ) for 45 s, and extension at 72°C for 45 s.A final step at 72°C for 7 min was used to complete primer extension.We included 120 replicate samples from independent PCR of a mix comprising the samples processed in this study and the samples that were provided from older studies (see above) to evaluate the reliability of sequencing.PCR products were checked on 1.5% agarose gel.Purification of PCR products was done with QIAquick or MinElute PCR purification kits (Qiagen) by pooling eight PCR products of similar amplification intensity per column (feasible due to the use of individual barcodes).Due to the large number of individuals used in this study, we pooled the samples into two libraries (A and B) according to the combination of individual barcodes and their equimolar concentrations.The concentration of pools of purified PCR products was quantified with Qubit 2.0 fluorometer (Life Technologies).The two libraries were sequenced with a 250 bp paired-end MiSeq protocol (Illumina) in one run by Fasteris (Switzerland).
MHC genotyping consists of applying two genotyping methods that we have previously confirmed to accurately provide MHC genotypes for barn owls based on allelic segregation within pedigrees (see Gaigher et al., 2016).We used the "Degree of Change" (Lighten et al., 2014) which identifies a coverage break point between the coverage of the last true allele and the next variant which would be an artefact.The variant with the highest breakpoint value is assumed to be the last true allele.We also applied the "Threshold method" (Galan et al., 2010) in which artefacts are separated from true variants based on the frequency of a variant within individuals.We defined a threshold at 6% and 9% for MHC-IIβ DAB1 and MHC-IIβ DAB2, respectively, since allele frequency distributions showed a substantial drop at these values suggestive of separation of true variants from artefacts.After both methods, we further inspected and manually removed variants with a chimeric composition comprised of parental true variants.Finally, we aligned the remaining variants and checked for the functionality of each gene with ClustalW in MegaX (Kumar et al., 2018).None showed evidence of non-functionality such as frame shift mutations or stop codons.

S2. Supertype clustering
Firstly, we checked for the presence of PSS on each MHC gene class separately using a maximum-likelihood analysis with the software CodeML implemented in PAML v4.8 (Yang, 2007).The input phylogenetic tree of each gene from MHC gene class and needed by CodeML was constructed with MrBayes under the GTR+I+G substitution model for MHC-Iα, while for both MHC-IIβ loci we used the JC+I+G substitution model, as inferred as the best-fit in jModelTest (Darriba et al., 2012;Guindon & Gascuel, 2003).We repeated the CodeML analysis with nine other randomly chosen trees as tree topology could bias the results (Gaigher et al., 2018).We compared two pairs of site models that contrast a neutral model of evolution vs. a model that allows sites' positive selection: M1a vs. M2a, and M7 vs. M8.M1 is the neutral model and assumes two classes of sites (conserved, ω = 0; neutral, ω = 1) in the protein.M2 adds a third class of site with ω as a free parameter, allowing for sites with ω > 1. Model M7 uses a β-distribution of sites between the intervals ω = 0 and ω = 1.M8 adds an extra class of sites to the M7 model, allowing for sites with ω > 1.The pairs of models are then compared using a likelihood-ratio test (LRT, given as 2*(lnLb -lnLa)) and compared to the χ 2 distribution with two degrees of freedom.The LTR revealed that for all genes, models that allowed for positive selection showed a better fit for both sets (Table S5).
Codons being under positive pressure were identified with the Bayes Empirical Bayes (BEB posterior probability > 95%) procedure (Yang, 2007).Secondly, all amino acids outside PSS were excluded (Minias et al., 2018) and then for each gene separately we characterized each PSS (see above) according to the physicochemical properties of each amino-acid site, based on five metric descriptors (z1, hydrophobicity;z2, steric bulk;z3, polarity;z4 and z5, electronic effects;Doytchinova and Flower 2005).This resulted in a matrix for each gene with the different MHC alleles in rows and z-values in columns (basically 1 row and 5 columns per allele).Using this matrix, alleles were clustered into supertypes by discriminant analysis of principal components (DAPC) with the adegent R package (Jombart, 2008).Firstly, we selected the number of clusters with the find.cluster()function with the number of clusters, k, varying from 1-30 for MHC-Iα and from 1-20 for each MHC-IIβ loci.We retained all principal components (PCs).We replicated this step 20 times (Migalska et al., 2019).After visual inspections of the plots of mean values of Bayesian Information Criterion (BIC) against the number of clusters, we selected the optimal number of clusters as the last consecutive increase in k to improve (i.e., reduce) BIC (i.e., typically referred as an "elbow" in the plots; Figure S1a).This procedure was replicated for 20 runs of clustering.With this approach, we did not find a stabilization point for MHC-IIβ and we proceeded with the "silhouette" method as an alternative method of clustering with the factoextra R package (Kassambara & Mundt, 2020) and the fviz_nbclust() function.The lack of stabilization point on BIC criterion was probably caused by a smaller number of alleles that are also more divergent in MHC-IIβ as compared to MHC-Iα.Although both methods use K-means algorithm to indicate the number of clusters present in the data, the difference between both methods is that the former uses intra-cluster distances while the latter uses a combination of inter and intra-cluster distances.While on the "elbow" method the best k is represented with the point where the BIC value is no longer improved, in the "silhouette" the best k is where the average silhouette value is the highest (Figures S2a and Figure S3a).Secondly, we then used the function DAPC to obtain clusters of the alleles.Then we chose the optimal amount of PCs to retain (this is, the PCs that explain the maximal variance on the dataset) with the function optim.a.score (n.rep=100) of the same package.Lastly, we repeated the DAPC using the optimal number of PCs and used the posterior probabilities to assign MHC-Iα alleles to the different groups (i.e.,supertypes;Figures S6b,S7b and S8b).S7 and S8 for best models included and model averaging results).Each plot shows the predicted slope for the father's (dashed blue) and mother's (solid red) MHC functional divergence and respective confidence intervals while holding other co-variables at their mean.

S3. Sample size
Presence/absence of supertypes Table S9.Summary of the models testing the influence of the presence/absence of specific supertypes on clutch size.Ages and laying dates were fixed in all models (see methods).Significance values (P value) were adjusted according to Benjamini-Hochberg procedure (Yoav & Yosef, 1995) to correct for multiple testing (P adj.).Fa: father; mo: mother; LD: laying date.S10 and S11 for best models included and model averaging results).Each plot shows the predicted slope for the father's (dashed blue) and mother's (solid red) MHC functional divergence and respective confidence intervals while holding other co-variables at their mean.The MHC functional divergence of either parent did not show evidence of explaining the observed fledging success and despite the apparent steep slopes, the wide confidence intervals convey the uncertainty of MHC divergence as an explaining factor of fledging success (see Discussion section).

MHC-Iα
Table S13.Post-hoc Tukey pairwise comparisons of the fledging success between all couples regarding the presence or absence of MHC-IIβ DAB1 supertype 5 in each parent."0": absence, "1": presence; "♂": father, "♀": mother.Tukey adjustment for multiple comparisons was applied and tests were performed on the log odds ratio scale.Results are averaged over the other covariates.S12 for simplified summary tables of models.Each plot shows the estimated marginal means and CI for each combination of parents regarding the presence of each supertype.Supertype 5 of MHC-IIβ DAB1 (third panel on B) shows a significant decrease in the fledging success when both parents carry it (Table S13).S14 and S15 for best models included and model averaging results).The right-hand of the plots shows the predicted slopes for the MHC functional divergence and their CIs for the genetic parents and the left-hand for the social parents.The MHC-Iα functional divergence of the genetic father shows evidence of explaining the observed fledging success even though his chicks are raised in another nest (panel A; Table S15).The genetic mother's MHC-Iα functional divergence is not shown since it is not included in the top models (Table S15).

S7. Genetic vs. social parents -fledging success
Table S17.Post-hoc Tukey pairwise comparisons of the fledging success between all genetic couples regarding the presence or absence of MHC-IIβ DAB1 supertype 2 in each parent."0": absence, "1": presence; "G♂": genetic father, "G♀": genetic female.Tukey adjustment for multiple comparisons was applied and tests were performed on the log odds ratio.Results are averaged over the other covariates.Refer to the methods section in the main text for modelling structure; and Table S16 for simplified summary tables of models.Each plot shows the estimated marginal means and CI for each combination of genetic parents regarding the presence of each supertype.Supertype 2 of MHC-IIβ DAB1 (second panel of B) shows a significant decrease in the fledging success when only the mother carries this supertype, as compared to when neither parent carries it (Table S17).The supertype 4 of MHC-Iα cannot be modelled with interaction due to low sample size, thus only the main effects are modelled and plotted.

Figure S1 .Figure S2 .Figure S3 .
Figure S1.Clustering of MHC-Iα alleles by Discriminant Analysis of Principal Components (DAPC), based on the physicochemical properties of translated amino acids inferred to comprise the Peptide Binding Region (PBR).a) Number of the most appropriate number of clusters given by the lowest BIC value after which it increases, in this case, it is K = 9. b) Clustering visualization retaining the first and second axes of the Discriminant analysis (DA) which are the ones retaining most of the variation observed (DA eigenvalues, bottom left).

Figure S4 .
Figure S4.Clutch size as a function of the functional divergence between the alleles within A) the MHC-Iα; B) the MHC-IIβ DAB1; C) the MHC-IIβ DAB2.The fitted regression line is obtained from the estimated marginal means of the averaged GLMM model (see methods in the main text for modelling structure, and TableS7and S8 for best models included and model averaging results).Each plot shows the predicted slope for the father's (dashed blue) and mother's (solid red) MHC functional divergence and respective confidence intervals while holding other co-variables at their mean.

Figure S5 .
Figure S5.Estimated marginal means from GLMMs analysing the effect of the presence/absence of specific supertypes of A) MHC-Iα, B) MHC-IIβ DAB1; and C) MHC-IIβ DAB2 on the clutch size.Refer to the methods section in the main text for modelling structure; and TableS9for simplified summary tables of models.Each plot shows the estimated marginal means and CI for each combination of parents regarding the presence of each supertype.

Figure S6 .
Figure S6.Fledging success as a function of the functional divergence between the alleles within A) the MHC-Iα; B) the MHC-IIβ DAB1; C) the MHC-IIβ DAB2.The fitted regression line is obtained from the estimated marginal means of the averaged GLMM model using the whole dataset (see methods in the main text for modelling structure, and TableS10and S11 for best models included and model averaging results).Each plot shows the predicted slope for the father's (dashed blue) and mother's (solid red) MHC functional divergence and respective confidence intervals while holding other co-variables at their mean.The MHC functional divergence of either parent did not show evidence of explaining the observed fledging success and despite the apparent steep slopes, the wide confidence intervals convey the uncertainty of MHC divergence as an explaining factor of fledging success (see Discussion section).

Figure S7 .
Figure S7.Estimated marginal means from GLMMs analysing the effect of the presence/absence of specific supertypes of A) MHC-Iα, B) MHC-IIβ DAB1; and C) MHC-IIβ DAB2 on fledging success.Refer to the methods section in the main text for modelling structure; and TableS12for simplified summary tables of models.Each plot shows the estimated marginal means and CI for each combination of parents regarding the presence of each supertype.Supertype 5 of MHC-IIβ DAB1 (third panel on B) shows a significant decrease in the fledging success when both parents carry it (TableS13).

Figure S8 .
Figure S8.Fledging success as a function of the functional divergence between the alleles within A) and B) the MHC-Iα; C) and D) the MHC-IIβ DAB1; and E) and F) the MHC-IIβ DAB2.The fitted regression lines are obtained from the estimated marginal means of the averaged GLMM model using the dataset subject to cross-fostering (see methods in the main text for modelling structure, and TableS14and S15 for best models included and model averaging results).The right-hand of the plots shows the predicted slopes for the MHC functional divergence and their CIs for the genetic parents and the left-hand for the social parents.The MHC-Iα functional divergence of the genetic father shows evidence of explaining the observed fledging success even though his chicks are raised in another nest (panel A; TableS15).The genetic mother's MHC-Iα functional divergence is not shown since it is not included in the top models (TableS15).

Figure S9 .
Figure S9.Estimated marginal means from GLMMs analysing the effect of the presence/absence of specific supertypes of A) MHC-Iα, B) MHC-IIβ DAB1; and C) MHC-IIβ DAB2 on the genetic parents on the fledging success.Refer to the methods section in the main text for modelling structure; and TableS16for simplified summary tables of models.Each plot shows the estimated marginal means and CI for each combination of genetic parents regarding the presence of each supertype.Supertype 2 of MHC-IIβ DAB1 (second panel of B) shows a significant decrease in the fledging success when only the mother carries this supertype, as compared to when neither parent carries it (TableS17).The supertype 4 of MHC-Iα cannot be modelled with interaction due to low sample size, thus only the main effects are modelled and plotted.

Figure S10 .
Figure S10.Estimated marginal means from GLMMs analysing the effect of the presence/absence of specific supertypes of A) MHC-Iα, B) MHC-IIβ DAB1; and C) MHC-IIβ DAB2 on the social parents on the fledging success.Refer to the methods section in the main text for modelling structure; and Table S16 for simplified summary tables of models.Each plot shows the estimated marginal means and CI for each combination of genetic parents regarding the presence of each supertype.Supertype 4 of MHC-Iα and supertype 1 of MHC-IIβ DAB1 cannot be modelled with interactions due to low sample size, thus only the main effects are modelled and plotted.

Table S1 .
Sample size of breeding attempts used per year for clutch size and fledging success analysis of each gene.Nb. obs: number of breeding attempts; ♀: number of different females; ♂: number of different males.Total values correspond to the total number of breeding attempts analysed and the total number of different females and males included.

Table S2 .
Sample size of different breeding attempts used per year for fledging success analysis of each gene on the cross-fostered dataset.Nb. obs: number of breeding attempts; ♀: number of different females; ♂: number of different males.Total values correspond to the total number of breeding attempts analysed and the total number of different females and males included.

Table S4 .
Frequency of MHC-IIβ DAB1 and MHC-IIβ DAB2 alleles in the barn owl population of western Switzerland and their GenBank accession numbers.New alleles found in this population are shown in bold.

Table S5 .
Summary of the analysis of positive selection.See the supplementary material and methods section for more details.

Table S7 .
Models with AICc ≤ 2 explaining the relationship between clutch size and functional divergence of the mother and the father at MHC-Iα, MHC-IIβ DAB1 and MHC-IIβ DAB2 genes.The age of the individuals, and first and second-order levels of laying date were fixed in all models (including the null model) and were omitted for simplification (see section 2.6).df: degrees of freedom; AICc: Akaike information criterion with correction for small sample sizes; ΔAICc: difference in AICc to the top model with the lowest AICc; w: Akaike weight for each candidate model, as the probability of each model being the best model in the set.ER: a measure of how much better one model explains the data than the next model.

Table S8 .
Summary of GLMMs analyses that tested the influence of functional divergence on the clutch size of barn owls.The standardized estimates are unconditionally averaged from the models within the top two units of AICc model ranking (TableS7).Estimates are standardized in two SE and presented with 95% Confidence intervals (lower and upper CI).The response variable (clutch size) has errors under the Poisson distribution with a logLink function.

Table S14 .
Models with AICc<2 for the relationship between fledging success and functional divergence of genetic and social parents at MHC-Iα, MHC-IIβ DAB1, and MHC-IIβ DAB2.The age of the individuals, and first and second-order levels of laying date were fixed in all models (including the null model) and were omitted for simplification (see section 2.6).df: degrees of freedom; AICc: Akaike information criterion with correction for small sample sizes; ΔAICc: difference in AICc to the top model with the lowest AICc; w: Akaike weight for each candidate model, as the probability of each model being the best model.ER: a measure of how much better one model explains the data than the next model.GFa: genetic father; GMo: genetic mother; SFa: social father; SMo: social mother.